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Abstract: We extend the massless dipole formalism of Catani and Seymour, as well as its 
massive version as developed by Catani, Dittmaier, Seymour and Trocsanyi, to arbitrary 
helicity eigenstates of the external partons. We modify the real radiation subtraction terms 
only, the primary aim being an improved efficiency of the numerical Monte Carlo integration 
of this contribution as part of a complete next-to-leading order calculation. In consequence, 
our extension is only applicable to unpolarized scattering. Upon summation over the 
helicities of the emitter pairs, our formulae trivially reduce to their original form. We 
implement our extension within the framework of Helac-Phegas, and give some examples 
of results pertinent to recent studies of backgrounds for the LHC. The code is publicly 
available. Since the integrated dipole contributions do not require any modifications, we 
do not discuss them, but they are implemented in the software. 
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1. Introduction 

At present, the necessity of next-to-leading order (NLO) calculations of QCD backgrounds 
for the Large Hadron Collider (LHC) is unquestionable. Much effort has been put into this 
problem, but until only recently, it seemed that the task is so huge that theory will stay 
behind the needs of experimentalists for quite some time. Whereas impressive calculations 
have been done with traditional methods based on Feynman diagrams [1-21], it is the new 
unitarity based methods [22-28] that provide some hope for accelerated progress. By now, 
there are three major groups with advanced software for virtual corrections [29-41], closely 
trailed by independent efforts [42]. Moreover, a full automate based on traditional methods 
is being built [43-46]. In any case, first successes have been recorded in [47-52] and look 
very promising. 

Any NLO calculation consists of two parts, which are separately infrared (soft/collinear) 
divergent, the virtual corrections and real radiation. In order to allow for Monte Carlo 
simulations, two general classes of approaches have been devised, namely phase space slic- 
ing [53-58] and subtraction [59-65]. Currently, the latter approach seems to have proven 
its superiority. Its most accepted version, dipole subtraction, has been presented in the 
massless case in [62] and later for massive partons in [65] (see also [66]). There are a few 
completely automated implementations of this method, which have been presented in the 
literature [67-70], but only one is available for public use [70]. Moreover, it has often been 
criticized that the large number of terms in the dipole subtraction formalism is a practical 



- 1 - 



problem in realistic calculations, since it leads to a high computational costs. With the 
present publication we want to remedy both problems. 

The basic idea is to allow for the same optimizations as those used in leading order 
simulations. The most important of these, besides phase space optimization, is the replace- 
ment of exact summation over external state polarizations by a probabilistic approach. In 
order not to ruin the Monte Carlo error estimates, the treatment must be consistent be- 
tween the real emission contribution and the appropriate subtraction. This is the main 
subject of the present publication. Note that partial extensions of the dipole subtraction 
formalism to polarized states have been discussed before in [71,72], but they were only con- 
cerned with scattering of polarized fermions or photons. The formulae from these studies 
are consistent with ours, when summed over the polarization of the gluon (or photon as 
discussed by these authors). Besides giving the appropriate formulae, we implement them 
within the framework of a fully fledged Monte Carlo generator, Helac-Phegas [73-75], 
which has, on its own, already been extensively used and tested in phenomenological stud- 
ies (see for example [76-79]). We demonstrate the potential of the software by performing 
some realistic simulations, in particular of the process gg tibbg^. We also argue that the 
inclusion of the subtraction terms does not increase the total evaluation time per phase 
space point by large factors. In fact, the actual source of the substantially longer eval- 
uation times in comparison to a leading order evaluation is the more complicated phase 
space, requiring orders of magnitude more accepted points to reach the same accuracy. 
In consequence, one can either further improve the evaluation time per phase space point 
by using color sampling for example (trivial in our approach), or concentrate on a better 
description of the phase space. We leave both tasks to future studies. 

The paper is organized as follows. In the next section, we discuss the role and treat- 
ment of polarization in Monte Carlo generators, which leads us to the motivation for the 
present study. Subsequently, we describe the behavior of cross sections in soft and collinear 
limits, when polarized partons are present. Section |^ contains our main results, namely 
the dipole subtraction formulae for polarized external partons. Section |5|, on the other 
hand, presents a few details of our implementation within Helac-Phegas, as well as some 
realistic simulation results. We conclude in Section 

2. Polarization in Monte Carlo simulations 

Most practical problems, which are solved with Monte Carlo simulations involve unpolar- 
ized particles. In consequence, it is necessary to sum/average over the spin of the incoming 
and outgoing states. This increases the computational complexity of a calculation by a 
factor, which can, in principle, amount to 2"'2 3"^, where n2 and are the numbers of 
particles with 2 and 3 polarization states respectively. There are usually some symmetries, 
like the chiral symmetry in the massless fermion approximation, or supersymmetry in the 
pure gluon case, which reduce the number of degrees of freedom. Moreover, with modern 
recursive methods for tree level matrix element evaluation it is possible to reuse parts of 

^While this publication was being prepared, a complete study of the full hadronic process has been 
published [21]. 
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a result to speed up the summation. Nevertheless, deterministic, exact methods have an 
inherent slow down factor, which cannot be completely removed. 

Since the phase space integration is already done with probabilistic methods, and 
the polarization sum is not coherent, it is clearly desirable to replace this sum by some 
kind of random sampling as well. The approach, which is most often used is to sample over 
helicity. The main disadvantage here is that different helicity configurations contribute very 
differently to the final result. In fact, several orders of magnitude between contributions 
are usually observed. If we use a fiat distribution to pick a helicity configuration, the final 
variance of a result will be increased by a substantial factor, not to mention the increase in 
the number of generated phase space points to obtain a given error estimate. The situation 
can be improved by using stratified sampling, which is usually also used for phase space 
optimization. One of the possible algorithms, which we implemented in our software, has 
been presented in [80] (more details will follow in Section |^ . 

Another approach, which is substantially easier to implement has been proposed in 
[81,82]. The idea is to replace summation over helicity by integration over a phase. For 
example, a gluon polarization state can be written as 

e^{k, (/>) = e'h^ik, +) + e-'^e^{k, -) , (2.1) 

where e^(A:,ib) are helicity eigenstates. The sum over the helicity of this gluon fulfills the 
following identity 

r-27r 



^|Ma|2 = ^ /"#|M^|2. (2.2) 

\ 



A 

Notice that the range of integration could have been reduced to [0, vr] with the same result. 
We keep 2tt as the upper bound in order to accommodate the third degree of freedom of 



massive gauge bosons, which is then added in Eq. 2.1 without any phase factor. Obviously, 
we do not expect large differences in the values of [Af,^!^, as function of (p, since for every 
value of (p we have both helicities contributing. In consequence, a fiat distribution in the 
Monte Carlo sampling should provide satisfactory results, which has been confirmed on 
specific examples in [81,82] and also in our studies. 

At this point, we should decide which approach to chose for dipole subtraction. As 
we will show in the next Section, helicity eigenstates provide particularly simple formulae 
for this problem, which are only minor modifications of the original formalism. Therefore, 
we trade the simplicity of the implementation of the Monte Carlo integration over general 
polarizations with a phase, for simple dipole subtraction formulae. Practice also shows that 
whereas helicity sampling is superior for a small number of helicity configurations, with 
many particles general polarizations with phases start to dominate. The reason is that 
at some point the algorithm of helicity sampling cannot find the optimal distribution. In 
our case, where the number of partons is relatively low, since we are always thinking of a 
next-to-leading order calculation, which is bound in complexity by our ability to compute 
virtual corrections, the large number of final states comes from electroweak decays into 
colorless states. We found it optimal to use a hybrid model, where partons have definite 
helicities and remaining particles have general polarizations with phases. 
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3. Soft and collinear limits for polarized partons 



In view of the considerations of the previous Section, we have two possibiUties to treat the 
polarization of partons. Indeed, we can either use arbitrary polarization vectors, or helicity 
eigenstates. Let us first show that, when using the latter the soft limit is particularly simple. 

It is well known that the exchange of a soft gluon between two partons (quarks or 
gluons) can be approximated by eikonal currents 

|M^+l|2~^(...|jA'.'^tja|...)^^ (3.1) 

where we have omitted irrelevant constants and the current J is given by 



J'^''^ = ^rf4:> (3-2) 



with k the momentum of the soft gluon and pi the momentum of a hard parton. The color 
charge operators are defined as in [62], i.e. for two given color space vectors |ai, . . . , am) 
and |6i, . . . , 6r),)) we have 

(ai, . . . ,ai, . . . , am\Tl\bi, . . . , 6j, . . . , hm) = ^aibi ■ ■ ■ ^aib, • • • '^flmfem ' (3-3) 

where = ifhac (color charge matrix in the adjoint representation) if parton i is a 
gluon, = (color charge matrix in the fundamental representation) if parton i is 
an outgoing quark, and T^^ = —t^^ if parton i is an outgoing anti-quark. The charges of 
in-going partons are defined by crossing. With this definition 

^T- = 0, (3.4) 



which means that no signs are needed in Eq. 3.2, and the current is both transverse and 
self-conjugate 

A;^ JM." = , J^'" t = J^^,'^ . (3.5) 

While it is clear that the eikonal approximation, and thus also the soft limit, is independent 
of the polarization of the hard partons, if the soft gluon is polarized on the other hand, we 
have 

j/.,atja _^ -J^^'-\r^\^{k,\K{kA) . (3.6) 

Crucially, for helicity eigenstates 

6*(A;,+) = e*'/'6^(A;,-), (3.7) 

where ■0 is some phase, which can be freely chosen. With this relation and the properties 
of the eikonal current, it is easy to show that 

J,,a t J-,a^^(^^ (^^ ^) ^ j^.,a t J-.a^^(^^ (^^ _) ^ ^_j,,a ^ ja (3 3) 

Thus, we have shown that the soft limit is independent of the helicity of the soft gluon. 
This is important, because this means that we can introduce helicity into the unpolarized 
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dipole subtraction formulae without explicit reference to the polarization vectors of the soft 
gluons. Unfortunately, for general polarization vectors, such as those defined in Section |2|, 



there is no relation of the type of Eq. 3.7. In consequence, if we wanted to use general 
polarization vectors, a less trivial modification of the dipoles would be needed. 

Let us now turn to the collinear limit. It is in principle possible to work directly with 
amplitudes, however due to the treatment of the soft limit, the original dipole subtraction 
formalism has been formulated for squares of matrix elements. Here we shall proceed 
similarly. 

We consider the quasi-collinear limit for massive partons and the corresponding true 
collinear limit in the massless case. For a pair {i,j} of outgoing partons, which become 
collinear, we assume the following momentum parameterization 



Pi 



(1 - - k'i 



z 2pn ' 

k\ + {l- zfmf- - m| 

\ — z 2pn 



(3.9) 
(3.10) 



where 



mt, pi- 



rn 



^ and 



n^ij, n is a light-like auxiliary vector, and k± is the 



transverse component orthogonal to both p and n, which parameterizes the collinear limit. 
The parton of mass rriij is the virtual particle, which splits into i and j. Its nature is 
uniquely determined in QCD. For example, if z is a quark and j is an anti-quark, then 
rriij = corresponds to a virtual gluon. As in [65], we define the limit by a uniform 
rescaling 

(3.11) 



k± ak± , rrii — > arrii , rrij arrij , 
with a ^ 0. The matrix element behaves in this limit as 



rrii 



am 



I] 1 



a+i(...,{pi,Ai},...,{pj,Aj} 
1 Svras 



m+l 



(3.12) 



2( _L ^m(---,{P,A'},...|Prj. .(2;,A;_L,{m},{A})|...,{p,A},...)m , 

o [Pi + Pj) m^j 



where P are generalized Altarelli-Parisi kernels (here in four dimensions), {m} is the set 
of masses, and {A} is the set of helicities, whereas ij is the emitter parton (we will call the 
original pair, the emitter pair). 

Whereas the unpolarized massive case of P has been presented in [65], the polarized 
massless case can be read to a large extent already from [83]. Here we present the formulae, 
which contain all the information 



{^'\PQQiz, k±,mQ, Xq, Xg)\X) 



1 - z 



m 



Q 



'PQPg 



{5\\Q5xQ\g - ^Iaq^AqaJ - + ^'xqxS^xxq - ^xxg)'^ 



(3.13) 



'^PQPg 



5v 



X'X 
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{M\PgQ{z, k^,mQ, Xq, Xq)\X) 



(3.14) 



Syx6' 



'A 



2^AgAg ik±-e{p,X')r{k^-e{p,X)) 
{jPQ+Pqf 



+ 



XXq 



{X'\Pgg{z,k^,Xi,Xj)\X) = 



Ca{5. 



5A'A 



"^AAi 




(3.15) 

e(p,A)) 



+<^A'A [<^AA, - <^AA,] (1-2^)} ; 

where = ^ — ^XiXj- Notice that the contents of the square brackets vanish upon 

summation over the heUcities of the emitter partons. Because the soft hmit remains un- 
touched with hehcity eigenstates, the above formulae alone can be used to modify the 
dipolc subtraction terms. An important difference to previous studies is that we do not 
use open Lorentz indices, but contract directly with polarization vectors. This is allowed 
by the transversality of the k±_ vector (and its counterparts in the final dipole subtraction 
formulae) and of the matrix element. Moreover, we used the phase conventions for the po- 
larization vectors as given in [73], which are also consistent with those of MadGraph [84]. 
Any difference in conventions should be relatively easy to compensate. 

Notice finally, that the initial state splittings can be derived from the above formulae 
by crossing. 



4. Dipole subtraction with helicity eigenstates 

The dipole subtraction formalism has been described to great extent in [62] and [65]. We 
do not repeat this discussion and assume that the reader is familiar with the main concepts. 
On the other hand, we give all the necessary formulae for a complete implementation in a 
numerical program, i.e. not only the modified splitting kernels, but also the momentum 
remappings. 

The starting point of our exposition is the subtracted real radiation contribution to a 
next-to-leading order cross section 

y"d$^(iM^+ii2-p), (4.1) 

where the sum runs over polarizations and colors, and an average over the initial state 
colors and polarizations is understood together with a symmetry factor for the final states. 
We have also omitted the jet functions in the phase space integration d<I>, since they are 
irrelevant at the moment. In this somewhat schematic expression, the dipole contribution 
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can be decomposed as 



^ = E E + E E + E E + E E ' (^.2) 

{i,j}kj^i,j a a,i j^i a,i b^a 

where i,j,k denote final states, whereas a,b initial states. A pair of indices corresponds 
to the emitter pair and a single index specifies the spectator. The fact that the sum in 



Eq. 4.1 runs over the difference of the real emission matrix element squared and the dipole 
subtraction contribution is not a coincidence of course, since we want to have a cancellation 
of divergences for each helicity configuration. The formulae that we present in the following 
guarantee that a contribution corresponding to a given helicity configuration of the partons 
is finite in soft and collinear limits. In principle, we could also give the formulae in a form, 
in which the finiteness would extend to a given color flow, but we refrain from this in the 
present work. 



Clearly, in each of the dipoles in Eq. 4.2 all the polarizations, but those of the emitter 



pair, must be taken over from the matrix element M^+i. Since the sums run over all 
partons and our formulae are only valid for helicity eigenstates, we have to require that 
all partons be in helicity eigenstates. On the other hand, there is no restriction on the 
polarization states of the remaining particles. 

Finally, let us stress that since we are only interested in the subtraction for the real 
emission, we work exclusively in four space-time dimensions. 

4.1 Final state emitter and final state spectator 

The single dipole contribution is 

T^ijM- ■ ■ ' fe' ^i}' • • • ' {Pj'^j}^ • • • ' {Pk, -^fc}, • • •) = -7 7 2" ('^•3) 

[Pi + Pj) ^ij 

X E^- • • ' "L^'^J' '^'i' • • • ' fe' ^k}, ■ ■ ■ \—=2~^^ij,k\ {Pij, A}, . . . , {pk, Afe}, . . .) , 

A', A «i 

where the momentum remapping is given by {Q = Pi + pj + Pk) 



^KQMPr+Pj?.rnl) 



K = „ K - "t^Q' + ' Q' ' (4-4) 



Q^^ J 2g2 



9t, = Q^-ru, (4.5) 

with A the Kallen function 

A(x, y, z) = + + — 2xy — 2xz — 2yz . (4.6) 
There are three cases to consider 
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Q Q{Pii ^i) + giPji Aj)) with rrii = ruij = mg, mj = and 

(A'|VQ5,fe|A) = STrasCp (4.7) 



1 - Zi{l - Vij^k) 



m: 



Q 



PiPj 



^x'x ; 



g Q{Pi, \i) + QiPj, Aj), with mi = mj = mg, rriij = and 



Vij,k 



Q 



2 2r+^- ipi+p^f 



25[ 



XiXj 



+h'x 



~(m) u ~(m) u 



~(m) u ~(m) y 

4 Pi - Pj 



^'xiX^ [SxXiZi + 6xxjZj - - 



+ ('^aa.-5U,)- + (<^aa,-4a,;t 



(4.8) 



eJl(Pu,A')ei.(pjj,A) 



zjj {Pi+Pjy 



• 9 ^ giPi^ K) + g{P3^ Aj), with rrii = mj = rriij = and 
(A'lV^^.felA) = 87ra,CA 



(4.9) 



X is 



X'X 



^AA,: 



+ 



26' - Kz+z^/2 



1 - ^i(i - yij,k) 1 - Zj{i - Vij^k) Vij,k 



+ 



AiA^ 1 
Vij,k PiPj 



~(m) a ~(rn) a 

zl 'pl - z'- 'pr 



~{m) y ~(m) y 
Z\ 'Pi - Z) 'pj 



e;l(Pii,A')e^(pij,A) 



+^ [<5aa,(1 - 2zi) + <5aa,(1 - 25,)] 1 ; 

'^%i,k ) 



where S'^.^. = 1 — Sx^Xj- In ef'Ch case, the content of the last square bracket vanishes upon 
summation over hclicitics of the emitter pair. The remaining terms have exactly the same 
structure as in [65] and are only modified by delta's in helicity. The variables Zi, Zj and 
Uij^k are defined as follows 

(4.10) 
(4.11) 



Zi 1 Zn 



PiPk 



yij,k 



pipk + Pjpk 

PiPj 



PiPj +Pipk+Pjpk 

The gluon splitting requires additionally 



~(m) ~ 1/, N ~(m) - 1.. N 

= ^« - - '^^j,k) , Zj =Zj--[l- Vij^k) 



(4.12) 
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where Vij^k is the relative velocity between the emitter pair with momentum pi + pj and 
the spectator with momentum pk, and can be expressed as 



_ J[2^^l + (1 - f^j - ^] - - - 

with 

/in = mn/VO^ . (4.14) 

Similarly, the relative velocity between the emitter with momentum pij and the spectator 
with momentum pj. can be written as 



A(l,4,/.i) 

''j'^ = 1 _ »2 _ „2 ■ (4-15) 

Moreover, the terms proportional to the free parameter k require the introduction of 
where Uij^j is the relative velocity between pi + pj and 



(1 - i4 - fi] - 4Yy\,, - 4^2 ^1 

'''' ~ a-f^l-f^j- ^X}m,,k + M 



(4.17) 



At this point it is important to comment on the freedom in the modification of the original 
formulae of [65]. Since we only require that these be recovered upon summation over 
helicities, arbitrary factors can be introduced as long as they do not spoil the collinear 
and/or soft limits. We have used this freedom to make the terms, which vanish upon 
helicity summation, resemble those that do not. Thus, for example, we have introduced 
the coefficient Vij^k/vij^k ™- front of the last square bracket in Eq. and similarly 



in Eq. [4.9| . A similar freedom exists in the treatment of the parameter k. As it represents 
a soft /collinear safe modification, it is unconstrained in the singular limits, and we are free 
to distribute it among different polarizations as we please. Here, we have chosen an even 
distribution among the four possible combinations. 

4.2 Final-state emitter and initial-state spectator 

The single dipole contribution is 

Vfj{...,{pi,Xi},...,{p,,X,},...;{pa,Xa},...) = \_ ^ } (4.18) 



T ■ T 



. . , {pij,X'}, . . . ; {pa, Xa}, ■ ■ ■ I ° 2 " ' ' ' fej' -^1' ' ' ' ! fe' ^a}, ■ ■ ■) , 



v 

A', A *J 

where 

_ VaVi + VaVj - PiPj + \ {nifj . . . . „j 



PaPi + PaPj 



(4.19) 
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and the momentum remapping is given by 



(4.20) 



#.=pf+^-(l-a;,,-„K- 



There are three cases to consider 



(4.21) 



• Q ^ Q{Pi, Xi) + giPj, Aj), with rrii = ruij = mg, rrij = and 



{SxXiSxiX^ - S'xx/x^xJt + ^'x.xASxXi - Sxx,)zi 



(4.22) 



^ Xi) + Q(pj, Aj), with nii = nij = mg, rriij = and 



(A'|V«q|A) = STTftsTR 



2(5', 



~ u ~ u 



^XiXj (^SxxA + SxxjZj - 



(4.23) 

[ziPi - ZjPj] el{pij, X')e,,{pij, A) 



m 



Q 



ZjJ [Pi+Pjf 



9 g{Pi-, Xi) + giPj, Xj), with rrii = mj = rriij = and 



X I Sx' 



(4.24) 



A'A 



hxi 



+ 



<^AA, 



1 •^i "l~ (1 Xij^a) 1 ~l~ (1 ■^ij,o, 



^^XiXj 



Pipj 



ZiPi - ZjP^- 



~Zip\ - Zjp'< 



^l{Pij,>'')^APij^>') 



+Sx'X [SxXi (1 - 2zi) + Sxx, (1 - 2zj)] } ; 



where 



PaPi 



PaPj 



(4.25) 



PaPi + PaPj PaPi + PaPj 

As in the previous case, terms in the last square bracket vanish upon summation over 
emitter pair helicities. 
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4.3 Initial-state emitter and final-state spectator 

The single dipole contribution is 

Vf{...,{pi,Xi},...,{pj,Xj},...-{pa,Xa},...) ^ ^ 



(4.26) 



^(...,{pj,Aj},...;{pai,A'}, 



T T 



^i^Vf|...,fe,A,},...;K,,A},...), 



A',A 



where 



PaPi + PaPj - PiPj 
PaPi + PaPj 



and the momentum remapping is given by 

P^ = pf + P^-{i-x,j,M, 



(4.27) 

(4.28) 
(4.29) 

As implicitly assumed above, we require both emitter masses to vanish, i.e. rrii = rrij = 
rriij = 0. There are four cases to consider 

• QiPa, Ao) g{Pi, Ai) + q, with 



Pai ~ ^ij,a.Pa 



(A'lVf |A) = 87ra,CF 



(4.30) 



X 



• ffCPa, Aa) q{pi,Xi) + q, with 
(A'lVf |A) = Stt^sTr 



(4.31) 



X {<^A„A,(1 - 2a;y- „(1 - Xij,a)) + [1 - 2(5A„Aj4,a} <^A'A<^AA, i 

• qiPa, Aa) ^ <7(Pi, Ai) + g, with 



(A'lVf |A) = SyrasCp 



{'^A'A^AAa^A^Aia^ij.a 

l-Xij^aUiC^-Ui) fpf 

1 — 1 /I i_ 

^ij,a PiPj \ -I- / \ 

+^A'A(^AaAi[2(^AAa - 1]} ; 



(4.32) 



e//(Poi,A')e*(pai, A) 



• 9{Pa, Aa) ^ ffCPi, Ai) + g, with 

'^AA„ 



(A'lVf |A) = Sna^CA 



X 5. 



'A' A 



1 Xij^a 



+ <^AAi (~1 + ^jj,a(l — ^^ij^a)) 



+<^AaAi 



1 - Xjj^a -», (! - Uj) f pf 



•^ij,a PiPj 



Pi 



Ui 1 — Ui I \Uj 1 — u. 



Pi 



(4.33) 



(■n{Pai,>^)el{pai,>) 
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where 



Ui = ; . (4.34) 

PiPa+PjPa 



4.4 Initial-state emitter and initial-state spectator 

The single dipole contribution is 

V-''\. . . , {pi, Xi}, {pa, Xa}, {Pb, h}) = -w^ — (4.35) 

X ^(~; {p,i,X!}, Afe}|^^^-^V-'^|~; A}, {p,, Afc}) , 

A', A 

where 

PaPb - PiPa - PiPb o«x 

Xi,ab = , (4.36) 

PaPb 

and the momentum remapping is given by 

K^ = ^^,abP^a, (4-37) 

~ 2kj-{K + K) ^ ~ ^ 2J^_K 

where the index j runs over aU final states and the momenta and are defined by 

K^^=p'a+P',-Pt , 

We require again both emitter masses to vanish, i.e. rrii = rrij = rriij = 0. There are four 
cases to consider 

• QiPa, K) g{Pi, K) + Q, with 

{X'\yi9'''\X) = STrasCp (4.40) 

X 



{t-TZ '^A„Ai(l+a;i,a6) MA'A-^AAa ; 

^ •'^i,ab ) 



• g{Pa, Aa) q{Pi, Ai) + q, with 

(A'|V^«'^|A) = Sna^Tn (4.41) 

X {'^A„A,(1 - 2Xi,„6(l - Xi^ab)) + [1 - 25A„Aj4aJ ^A'A<^AA, i 

• q{Pa, Aa) q{Pi, Ai) + g, with 

(A'|V««'^|A) = STrasCp (4.42) 

X {<^A'A^AA.'^A.Aia:i,a6 

1 ~ ^i,ab Pii ' Pb ( ^ PiPa fi\ ( „v PiPa ^i/ 



+Kxr—^ ^i^L^^ - ^-^pi \\p\- ^pi 

Xi,ab Pi -Pa Pi -Pb V PbPa J \ PbPa 

xe^(^aj,A')e*(pai,A) +(5A'A<5AaAi[2(5AAa - 1]} ; 
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• 9{Pa, Aa) 9{Pi, Ai) + g, with 



X ISx'x - 




(4.43) 




SxXa)'^Xi^ab] 



5. Implementation and example results 

As part of the present study, we have implemented the complete formalism given in the 
previous section within the framework of Helac-Phegas. The software can be obtained 
from the Helac-Phegas web page [85]. The main features can be summarized as follows 

1. Arbitrary processes 

We use the matrix element generator of Helac-Phegas, which means that any 
process, which can be calculated by this generator is also accessible to the dipole 
subtraction software. The only limitation is given by the models implemented. Cur- 
rently, the full Standard Model, both electroweak and QCD, is available. Let us 
stress that the subtraction is only applied to partons, and the software is therefore 
devised for NLO QCD corrections. The observables are specified by user defined jet 
functions and histogramming routines. The present implementation contains a built 
in kx jet algorithm. 

2. Massive and massless external states 

We have implemented the formulae as given in the previous section. This allows us 
to treat massive and massless partons on the same footing. As noted in the previous 
section, we assume the initial state partons to be massless. 

3. Helicity sampling for partons 

The integration over the phase space is done in three stages. At first, the phase space 
is optimized using multi-channel methods. In this phase, the user can either evaluate 
the subtracted matrix elements with full summation, or use a fiat helicity Monte 
Carlo, where all helicity configurations occur with the same frequency. Since this part 
is only used to determine the weights of the different channels we recommend to run in 
the latter mode, which is fast and gives exactly the same results as full summation. 
In a second stage, the phase space measure is fixed and the helicity configuration 
sampling weights are determined by evaluating the subtracted matrix element for 
all helicity configurations. The optimization consists of a standard minimization of 
the variance [80], but is as slow as full summation (in fact, this is full summation). 
We recommend to run this phase on a few hundreds up to a few thousand points. 
The optimum depends strongly on the number of helicity configurations. Finally, in 
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the last stage, the matrix element is evaluated with helicity sampling, which updates 
itself at user specified intervals in the number of accepted events. These intervals 
should be long enough, to allow the channels which have a low weight to accumulate 
enough events. We recommend a number of the order of ten thousand. All the 
parameters can be changed by the user. Obviously, one can make a complete run 
using full helicity summation of flat Monte Carlo. As with any simulation, some 
experimentation is needed to obtain best results. 

4. Random polarizations for non-partons 

Non-partons can be treated in two different ways as far as their polarization is con- 
cerned. The user can either require random phases as described in Section |2|, or 
treat non-partons on the same footing as partons. We recommend to use random 
polarizations. 

5. Restrictions on the subtraction phase space 

We have implemented a restriction on the phase space of the subtraction as proposed 
in [6]. This amounts to only including dipole subtraction terms, which satisfy the 
following criteria 

(a) final-final dipoles 

yij,k < a^ax ; (5-1) 

(b) final-initial dipoles 

(c) initial-final dipoles 



IF ( \ 

Ui < a^^^ ; (5-3) 

(d) initial-initial dipoles 

~ _ PaPi II , /(. A\ 

Vi — < ^max 1 y^-^) 

PaPb 

where the occurring variables have been defined in the previous section. The four 
c^max parameters can be freely chosen by the user. A value of 1 amounts to no 
restriction. Best results are obviously obtained with small Omax) since the matrix 
element for this setting behaves largely as the real emission alone and the phase 
space optimizer has been designed to emulate such behavior. Besides the phase 
space restriction, we have also included a technical cut. Using the above formulae 
one can evaluate a minimal value of the four parameters given on the left hand side 
for all dipoles. Let us denote it by amin- We reject a phase space point completely if 

"mm < acut ■ (5.5) 
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The value of Ocut should be specified by the user. Let us stress that the technical cut 
is necessary to avoid numerical instabilities in the cancellation between real emission 
and the dipoles. However, there are also cutoffs in the phase space generator, which 
are needed to avoid numerical instabilities in the generator itself. 

6. Phase space integration 

The integration is completely controlled by Phegas [74], which is a multi-channel 
phase space generator based on Feynman graphs. The only modification, which was 
necessary was the simultaneous treatment of both positive and negative weights. In 
fact, the user can specify whether only positive, negative or both weights should be 
included. This does not modify the histograms, however, which always contain the 
full result. 

We should mention that in the early stages of the development, we have adapted the 
kinematics remapping routine from MCFM [86]. Although, a minute part of the final code, 
it was important for debugging to be able to rely on the correctness of the kinematics. 

This software has been tested against MadDipole [87]. We would like, however, to 
point out that instead of comparing arbitrary processes, we concentrated on the individual 
parts of the code. In fact, the only sensitive code, which cannot be tested against the 
official version of Helac-Phegas is contained in the kinematics remapping, the color 
correlators and the dipole subtraction formulae themselves. In consequence, to ensure the 
correctness of the implementation it was sufficient to compare the latter for all independent 
cases. These are summarized in Tab. |l|. Of course, we could only compare the dipole 
subtraction terms after summation over helicities. The expressions for independent helicity 
configurations have been tested by checking for cancellation in the appropriate limits. 

Additionally, in Tab. |2| we have presented the measured time needed to evaluate the 
real emission matrix element and the subtraction terms. We note that the inclusion of 
the full set of subtraction terms slows down the computation by a factor of about three 
to four in most cases. However, since there is a restriction in the subtraction phase space, 
we expect (and indeed observe in practice) that the true additional cost of the subtraction 
does not exceed the cost of the real emission itself. Therefore, an improvement in the 
speed of the calculation will not be obtained by replacing the dipole formalism by another 
subtraction. In fact, a further substantial speed up can be obtained by turning to Monte 
Carlo summation over color configurations. What is the source of the high cost of the 
evaluation of the subtracted real radiation in comparison to leading order real radiation? 
It lies in the phase space integration. The integrand has now a much more complicated 
behavior, and therefore multi-channel optimization based on Feynman graphs does not lead 
to such a drastic improvement of the convergence. 

As a final demonstration of the power of our implementation, we have performed a 
simulation of the gg — ttbbg subprocess, which is part of the complete NLO calculation of 
pp -^■ ttbb+X. We have used the same set up as [18]. In particular, we have taken ^/s = 14 
TeV as center of mass energy and a top quark mass of mt = 172.6 GeV. The b quark has 
been kept massless. In order to obtain jets, we have used the fcy algorithm [88-90] with jet 
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So - massless emitter, So - massless spectator, £m - massive emitter, Sm - massive spectator, £i - initial 
state emitter, £p - final state emitter. Si - initial state spectator, Sf - final state spectator, / - check, ■ - 

does not occur. 





£o/So 


£o/Sm 


£m/So 


£m/Sm 


So/So 


£o/Sm 


£m/So 


£m/Sm 


£i/Si 










£f/Si 










n — ^ nn 


/ 


■ 


■ 


■ 


g^ gg 


/ 


■ 


■ 


■ 


g^m 


/ 


■ 


■ 


■ 


g ^qq 


/ 


■ 


/ 


■ 


q ^ gg 


/ 


■ 


■ 


■ 


q^qg 


/ 


■ 


/ 


■ 


q^gg 


/ 


■ 


■ 


■ 


q^gq 


/ 


■ 




■ 


£i/Sf 










£f/Sf 










g^gg 


/ 


/ 


■ 


■ 


g^ gg 


/ 


/ 


■ 


■ 


g^qq 


/ 


/ 


■ 


■ 


g^qq 


/ 


/ 


/ 


/ 


q^qg 


/ 


/ 


■ 


■ 


q^qg 


/ 


/ 


/ 


/ 


q^gq 


/ 


/ 


■ 


■ 


q^ gq 


/ 


/ 


/ 


/ 



Table 1: Independent dipole splitting formulae, which need to be tested in order to ensure the 
correctness of the code. In the splitting description, e.g. g gg, the left hand side particle always 
denotes the virtual state. 



Process 


Real Emission + Dipoles 


Real Emission 


Nr Of Dipoles 




[msec] 


[msec] 




gg ggg 


3.8 


1.0 


27 


gg gggg 


8.5 


2.6 


56 


gg -» ggggg 


300 


42 


100 


ud W+gggg 


9.3 


2.4 


56 


gg -» ttbbg 


12 


2.9 


55 



Table 2: The CPU time needed to evaluate the real emission matrix element together with all 
of the dipole subtraction terms per phase-space point (this corresponds to amax = numbers 
have been obtained on an Intel 2.53 GHz Core 2 Duo processor with the Intel Fortran compiler using 
the -fast option. 



-16- 




Figure 1: Distribution of the invariant mass mi^i of the bottom- anti-bottom pair (a), distribution 
in the transverse momentum pT^:^ of the bottom- anti-bottom pair (b) and distribution in the rapidity 
UijI of the bottom- anti-bottom pair (c) for pp{gg) ttbbg -\- X at the LHC for two different values 
of cxmax'- 0.01 (red solid curve) and 1 (blue dashed curve). 

recombination of partons with a pseudo-rapidity of \r]\ < 5 with AR < 0.8. Additionally, 
we required the b jets to satisfy pT,b > 20 GeV, \yb\ < 2.5. The phase space of the top 
quarks has not been restricted. The non-perturbative input of our computation has been 
given by the CTEQ6M PDFs [91,92]. Fig. |^ (a-c) contains the invariant mass, transverse 
momentum and rapidity distributions of the bb pair for Omax = 0.01 and Omax = 1 (all 
parameters set to the same value). The impact of the variation of amax is cleary visible. 
In fact, for both values, the shape resembles that of the leading order result, besides the 
normalization. On the other hand, for larger values of amax a relatively small negative dip 
for low invariant masses can be observed in the m^^j distribution. The latter is due to the 
strong enhacement of the dipole contributions for low m^j. 
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6. Conclusions 



We have presented an extension of the dipole subtraction formahsm to arbitrary hehcity 
eigenstates. Not only did we provide appropriate formulae, but also a public implementa- 
tion, which contains both the subtracted real radiation and integrated dipoles. The results 
of our tests show that with this software, calculations, which are of current phenomeno- 
logical interest, can be performed fully automatically within tractable time. In the nearest 
future, we will apply this technology to the processes from the "NLO wishlist" [93] using 
the automated 1-loop extension of Helac [41] for the virtual corrections. 
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